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In this paper we fill the gap in understanding the non-stationary resonance dynamics of the weakly 
coupled pendula model, having significant applications in numerous fields of physics such as super¬ 
conducting Josephson junctions, Bose-Einstein condensates, DNA, etc.. While common knowledge 
of the problem is based on two alternative limiting asymptotics, namely the quasi-linear approach 
and the approximation of independent pendula, we present a unified description in the framework 
of new concept of Limiting Phase Trajectories (LPT), without any restriction on the amplitudes of 
oscillation. As a result the conditions of intense energy exchange between the pendula and transition 
to energy localization are revealed in all possible diapason of initial conditions. By doing so, the 
roots and the domain of chaotic behavior are clarified as they are associated with this transition 
while simultaneously approaching the pendulum separatrix. The analytical findings are corrobo¬ 
rated by numerical simulations. By considering the simplest case of two weakly coupled pendula, 
we pave the ground for new opening possibilities of significant extensions in both fundamental and 
applied directions. 


The model of coupled pendula and some of its mod¬ 
ifications play a significant role in Mechanics pQ, Solid 
State Physics [2] including superconducting Josephson 
junctions EH3, Photonics, including Bose-Einstein con¬ 
densates [6], Biophysics, including DNA functioning [7]. 
The majority of the results in all these fields relate to 
stationary dynamics and are based on the fundamental 
stationary regimes, namely the Nonlinear Normal Modes 
(NNMs) in finite systems [8j |9] and solitons (breathers) 
in infinite models mm- As for non-stationary processes, 
NNMs can also be used for their description provided 
that the intermodal resonance is absent m- However, 
the non-stationary resonance dynamics of finite systems 
turns out to be much more complicated. Therefore only 
isolated and predominantly numerical results were ob¬ 
tained in this field mm- The recently developed con¬ 
cept of Limiting Phase Trajectories (LPTs) allowed for 
a systematic approach to description of non-stationary 
resonance regimes mans], including coupled pendula 
dynamics m- This concept introduces a fundamental 
non-stationary process of new type which corresponds to 
maximum possible energy exchange between the oscil¬ 
lators (in particular, pendula) or clusters of oscillators. 
In essence, the LPTs play in the non-stationary reso¬ 
nance dynamics of finite systems the role similar to that 
of the NNMs in the stationary theory and in the study 
of non-stationary yet non-resonant regimes. In terms of 
LPT, the transition from intense energy exchange be¬ 
tween some clusters of oscillators (coherence domains), 
in particular, weakly coupled pendula, to energy local¬ 
ization in the initially excited cluster (oscillator) can also 
be predicted DUES ED. 


However, all existing analytical results in the non- 
stationary resonance dynamics of finite-dimensional sys¬ 
tems relate to some asymptotic limits which are either 
quasi-linear systems or separated oscillators (pendula) 
HD EH- In this study we remove these restrictions and 
admit arbitrary oscillation amplitude in the framework of 
LPT concept. Assuming weak coupling, only the close¬ 
ness to inter-pendula resonance is required. It is shown 
that such an extension is crucial for revealing the nature 
of large amplitude dynamics and prediction of necessary 
conditions for the onset of chaotic regimes. The ana¬ 
lytical findings are confirmed by numerical simulations. 
We construct Poincare sections for the starting equations 
of motion to verify the obtained analytical results con¬ 
cerning with the revealed dynamical transitions and their 
connection with manifestation of chaotic behavior. 

For the sake of clarity we will discuss a mechanical 
interpretation of the problem in terms of two weakly 
coupled pendula undergoing planar motion as sketched 
in Fig.l. The corresponding dimensionless equations of 
motion can be written as 

dff n ■ 

+ sin q, + - q 3 _j) =0, j = 1,2 (1) 

where qj is the angular coordinate of the j-th pendu¬ 
lum, To = cjoL cjo = y/g/l is its linear natural fre¬ 
quency, g is the gravitational acceleration, 5 <C 1 and 
/3 is the coupling parameter. As known, these equations 
describe also a particular case of the two Josephson junc¬ 
tions (two-junction interferometer) |3] as well as of the 
Frenkel-Kontorova model, having numerous applications 
in solid state physics and photonics mm- 
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FIG. 1. Two identical pendula, with unit mass ra, weakly 
coupled by a torsional spring. 


Dealing with pendula oscillations under internal 1:1 
resonance conditions, we rewrite @ in the form 


cPq • 

-j^-+u 2 qj+e/3(qj-q 3 -j)+eij,(smqj-aj 2 qj) = 0 , j = 1,2 

( 2 ) 

where 0 < uj < 1 (the lower limit will be later on spec¬ 
ified) is the resonance oscillation frequency. Under in¬ 
ternal 1:1 resonance conditions the combination of two 
terms in the second brackets has to be small (we sup¬ 
pose of order- 5 ), and fi — s~ x is a book keeping param¬ 
eter. Thus, we assume the closeness to resonance but 
we don’t impose any restriction to the oscillations am¬ 
plitude and the ensuing resonance frequency. Passing 
to complex variables, given by ipj = + iujqj)e LUJT °; 

(p* = (^ — iuoqj)e~ luJTo , and substituting in (J2|, the two 
scale expansion pj = y^oCro, t\) + £^j,i(to,ti) + ..., in 
which t i = st o is slow time scale, then taking into ac¬ 
count that -£jr = + s-J^r and selecting the terms of 

order s° and e 1 , we get 


9< PJ ,o _ n 

<9 t 0 


difij, i difjfi 


dro dr\ 


+ ^ - y (</>?,o - l Pjfl e ~' 2zuTo - Vz-ifi - ft-j, 0 e- 2lMTo ) + H 


^., oe ^To _ 0 e-<u.T, 

2iuj 




— 2 iujto \ 


= 0, jm 1,2 


( 3 ) 


( 4 ) 


r 


The partial differential equations ([3| show that (pj^ = nary differential equations with respect to fast time scale 
(Pj,o(n) therefore equations Q can be considered as ordi- tq. While integrating them, the conditions of absence of 


dipjfl i/3 


dr\ uj 


-(<Pj,0 — 9?3—j,o) + p 




secular terms lead to the following equation 

a n / |^,o |2N " 


n= 1 


W 2 n +1 l g ) Vjfi + y ^-,0 


= 0, j — 1,2 (5) 


r 


in which a n = 1/(2 a n ) and a n is obtained by the follow- main asymptotic approximation in slow time scale with 

ing recurrence relations a n = b n a n - 1 , with ao = 1, and respect to (pj t o- Contrary to original system (1) which 

b n = 6 n _i + n, with 6 0 = 0. Equations © represent the possesses only one integral, system (J5| admits two inte¬ 
grals of motion: 


H = 


i/3 (iuj i 

~LO + ^ \T ~~ 2 w 


2 


3 = 1 


E l^'- 0 | 2 - -(^,0^2,0 + ^2,0^1,0) + E EF 1 ) 

— cj — z — 


2 oo 


n 


(l^o| 2 ) n+1 (6) 


j=l n=l 


(n + l)cj 2n+1 


^ = Em 2 

i=i 


( 7 ) 
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Taking into account the second integral 0, we set then, the first integral 0 , after series summation, takes 
^i,o = VNcosO(ri)e lSl ^ Tl \ </?2,o = VNsmO(ri)e lS2 ^ Tl>) ; the following form 


H— — cos A sin 20 H- 

ca 2 Nuj 


{A— 8<^ 2 +4 u; 2 [Jo(&2 cos#)+Jo(&2 sin#)]}, 


( 8 ) 


where A = 82 — £1 and k<i = y/~N/uo. Corresponding 


equations of motion can be written as follows 


0 = — — sin A 
uj 

A sin 2 0 = cos 20 cos A — -^= [ J\ (fc 2 cos 0 ) sin 0 — J\ (&2 sin 0 ) cos 0 ] , 


(9) 

( 10 ) 



FIG. 2. Evolution of the A -0 phase portrait for uj = 0.65 
and maximum angles qj,max — 37r/4 for decreasing parameter 
e; LPTs (red), separatrix (blue), a) Before first transition, 
£ sss 0.2; b) after first transition, e — 0.1; c) second transition, 
e — 0.0695; d) after second transition, e — 0.05. 


where Jq(), Ji() are Bessel functions of the first kind 
and the derivative with respect to t\ is considered. This 
system describes the slow dynamics of weakly coupled 
pendula for arbitrary initial conditions. Having obtained 
the integral 0 and equations of motions ([9])-([Tq]) we have 
the possibility to perform the analytical study of the 
considered problem. We begin by analyzing the phase 
portrait described by equation H= const.. It is conve¬ 
nient to present the evolution of the A -0 phase plane 
with changing the parameter £, characterising interpen¬ 


dulum coupling, for a given value of resonance frequency 
that determines also the maximum value Qj^max of pen¬ 
dulum oscillation angle qj. In Fig.2 the phase portrait 
corresponding to uj = 0.65 ( qj,max — 37r/4) and (3 = 1.0 
is shown. 

There exist two dynamical transitions, relating station¬ 
ary and highly nonstationary dynamics, that are mani¬ 
fested by changing the parameter s; the phase portraits 
shown in Fig.2 highlight the topological changes asso¬ 
ciated with these transitions. The stationary points in 
Fig.2 correspond to NNMs of the considered system; the 
closed phase trajectories surrounding them imply weak 
energy exchange between the pendula. The Limiting 
Phase Trajectory encircling all trajectories corresponds 
to full energy exchange between the pendula and the 
trajectories close to it describe intense energy exchange. 
The first dynamical transition occurs due to instability 
and bifurcation of the in-phase NNM (Fig.2b). It leads to 
appearance of two additional stable in-phase (yet asym¬ 
metric) NNMs, and the homoclinic separatrix encircles 
them. This is a local transformation of the phase portrait 
which strongly influences the stationary dynamics of the 
pendula. However, there is not yet any qualitative change 
in its highly non-stationary dynamics because complete 
energy exchange between the pendula can still be possi¬ 
ble. Only the second dynamic transition (Fig.2c), caused 
by global transformation of the phase portrait after co¬ 
alescence of LPT and homoclinic separatrix leads to a 
drastic change of non-stationary dynamics which can no 
longer give rise to full energy exchange between pendula 
and therefore is characterized by predominant energy lo¬ 
calization in the initially excited pendulum. Moreover, 
after this transition the homoclinic orbit turns into a het- 
eroclinic one (Fig.2d). 

The analytical conditions for both transitions are re¬ 
ported below and are confirmed by the numerical solution 
of equations 0. For the first transition prediction we re- 



















































4 


sort to the solution of equations © -flip]) in the vicinity 
of the stationary point A = 0, # = 7 t/4 (in-phase NNM) 
as A = Ai, 0 = 7r/4 + 6% Assuming that Ai and 0\ are 
small perturbations, the solution can be determined by 
means of the linearized version of equations 


Oi 

Ai 


-5a. 


c u 


K 

£Jo 

LJ 

UJ 

- J2 ( 
w V 

k 2 \ 

V2j 


M 

V2j 

Oi 



From the system (11)-(12) it is seen that 
in-phase NNM occurs when the coefficient 
to zero. The latter condition leads to the 
pression for the first transition 


instability of 
of 6 \ is equal 
following ex- 



(13) 

after which the phase portrait becomes qualitatively alike 
the one shown in Figure 2b. To reveal the condition of 
the second transition we derive the equation describing 
LPTs by considering that they possess the point 0 = 0 or 
6 = 7r/2, thus the corresponding Hamilton function can 
be written as Hlpt = H — C, in which C is the constant 
given by H for 0 = 0. From Hlpt and the equations 
of motion © -([10} the second order nonlinear differential 
equation for 0 , which is valid for LPTs only, is obtained. 
The corresponding (0, 0) phase portraits are presented in 
Fig.3b for several initial conditions. While the phase tra¬ 
jectories out of the separatrix correspond to full energy 
exchange between pendula, those inside the separatrix 
correspond to maximum possible energy exchange in the 
condition of energy localization in the excited pendulum. 
If the threshold of first transition is given by ( [13] ) (curve 
I in Fig.3a), the condition of second transition can be 
found by taking into account that its occurrence implies 
that LPT possesses the unstable stationary point A = 0, 
0 = 7r/4, therefore we get the sought second transition 


(curve II in Fig.3a), namely 


6 


2c j 2 

~W 


1 +Jo (fe) — 2 Jq 



(14) 


It is worth emphasizing that the described scenario is 
observed not only for small angles (quasi-linear case) 
but also for values of the angles close to 7 r. However, 
the threshold values of the parameters corresponding to 
both transitions change strongly with the resonance fre¬ 
quency, as shown in Fig.3a. The perfect agreement be¬ 
tween analytical prediction of curves /, II and their nu¬ 
merical counterparts /*,//* is observed for 0 < q < 37t/4 
(0 < cj < 0.65). The top horizontal axis labels refer to 
the amplitude q corresponding to the values of resonance 
frequency c 0 reported in the bottom horizontal axis. It 
can be seen that the lowest value of the latter resonance 
frequency is equal to c 0 = 0.2 (g maa; = 3.14). Moreover 
both boundaries shown in Figure 2a refer to parameter 
£ ranging from 0 to 0.2, in agreement with our initial 
assumption concerning with the smallness of the sum in 
the second bracket in equation Q. The quantitative dif¬ 
ference for larger angles (which reaches at most 10% and 
20% for curve I and //, respectively) can be reduced 
by considering next order approximation in the multiple 
scale expansion procedure. The structure of the phase 
plane depicted in Fig.3b allows to predict the temporal 
behavior of the angle variable 0. The trajectories situ¬ 
ated far from the separatrix correspond to almost straight 
lines. However, due to the restriction 0 ^ 0 ^ 7t/ 2, they 
become saw-tooth type functions. The analytical solu¬ 
tion of the problem in terms of non-smooth functions can 
be obtained after change of temporal variable through the 
procedure proposed in [22], 23 and used for the study of 
non-stationary resonance processes in P3JEE5JEE6]. As for 
phase trajectories located inside the separatrix, they cor¬ 
respond to localized LPTs and can be easily found after 
linearization of the second order equation for LPTs in 
the vicinity of 0 = 6 = 0. By substituting 6 = 0 in (J8| 
we get H = C , where C is the constant 


Li(Au 2 J 0 (k 2 ) + N -4co 2 ) 

2Nu 

Thus, Hlpt = H — C and the corresponding expression 
for cos A reads 


cos A = 


liuo 2 esc (0) sec (0) (—Jo (sin^)/^) — Jo (cos^)/^) + Jo (^ 2 ) + 1) 

7V/3 


(15) 


Expression (15) can be substituted into the first deriva- tive of given by 


6 = — — cos A A 


UJ 


(16) 
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FIG. 3. Evolution of the dynamic transitions in the (u — e) parameter space. The upper horizontal axis shows the maximum 
angle q for corresponding uj in the lower axis, a) first and second dynamic transitions, analytical prediction (I,II) and numerical 
observation (I*,II*); b) LPTs phase plane; c) temporal behavior of 0 corresponding to the trajectories shown in b). 


and A is obtained from 0. leading to the sought second 


order differential equation 


423/2 

0 = ^' esc 2 20 [-Jo (sin 0k 2 ) ~ Jo (cos 0k 2 ) + Jo (fe) + 1] 

— VN cosOJi (k 2 y/ujsm0) + VN sinOJi (k 2 ^/ui cos 0) + 2^/uoJo (k 2 ) cot 2 6 
—2cot(2$) Jo (k 2 cos 6) — 2\/Zj cot 20 Jo {k 2 sin 0) + 2cot 20] (17) 


Equation allows to construct the LPTs phase plane 
shown in Figure 3b. 

It is worth emphasizing that Fig.3a clarifies the relation 
among the obtained results and conventional approxima¬ 
tions used for description of coupled pendula dynamics: 
the quasi-linear approaches can be applied only in the 
right part of the parametric plane (0 < q < 7t/4) whereas 
the independent pendula approximation holds only for 
the bottom part of the parametric plane (0 < 5 < 0.03). 
On the contrary, the proposed approach based on reso¬ 
nance asymptotic turns out to be valid for description of 
regular motion in all parametric plane (u-e). 

So, we have found the dependence on the resonance 
frequency (or initial angle) of the thresholds correspond¬ 
ing to both dynamical transitions. These relations allow 
to single out in parametric space the domains of regular 
motion of the pendula. It has been shown that all regu¬ 
lar motions revealed in main asymptotic approximation 
are observed also by integrating the starting equations of 
motion 0 . Moreover, the analytical predictions of both 
dynamic transitions are well confirmed as well. However, 
it must be underlined that, contrary to the asymptotic 
approximation, the initial system 0 is not integrable. 
Therefore it is of interest to clarify the onset of chaotic 
behavior and the role played by LPTs in the general 


behavior of the pendula. Towards this goal, Poincare 
sections constructed on the basis of the starting equa¬ 
tion of motions 0 are reported in this section. In Fig.4 
and Fig.5 Poincare sections are shown referring to max¬ 
imum angles q^max — 37 t /4 and q^max — 9tt/ 10, respec¬ 
tively. The four sections correspond to different dynamic 
regimes (see Fig.2) for decreasing values of e, according 
to the points highlighted in Fig.3a, Fig.4a and Fig.5a re¬ 
fer to the dynamics before the first transition where the 
LPTs (red curve) encircling the in-phase NNM are also 
depicted. Fig.4b and Fig.5b refer to the case in-between 
the two transitions, where the new stationary states born 
as a result of instability of in-phase NNM can be seen; 
moreover the associated homoclinic separatrix encircles 
the corresponding stationary points. Fig.4c and Fig.5c 
reflect the conditions at second transition, where LPT 
becomes separatrix; as indicated in Fig.3a, manifesta¬ 
tion of chaotic behavior can be observed in the vicin¬ 
ity of second dynamic transition for large enough angles 
( qj,max — 97 t / 10 , Fig.5c). In Fig.4d and Fig.5d the lo¬ 
calized LPTs as well as the heteroclinic separatrix (blue 
curve) are well seen; it is worth noticing that, for small 
enough £, the motion remains regular in all phase-space 
(see also Fig.3a). 

Examples of the pendula oscillations temporal evolu- 
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tion corresponding to the parameters considered Figures 
4,5 are reported in the following Figures 6,7. The pen- 
dula response reported is obtained by direct numerical 
integration of the initial dimensionless equations of mo¬ 
tion §. More specifically, Figure 6 refers to the case 
uj = 0.65 and the maximum angles qj are 3/47T for which 
only regular behavior was observed (see Fig. 3a). Differ¬ 
ently, in Figure 7 the case uj = 0.482 and the maximum 
angles qj are 9/107T is considered. In agreement with the 
findings shown in Fig 3a, chaotic behavior can be ob¬ 
served in Figure 4c in the vicinity of second tansition. 

While only regular motion is observed for angles less 
than 135°, for larger angles signs of chaotization appear 
in the vicinity of the second dynamic transition triggered 
by heteroclinic chaos, as anticipated in the discussion of 
Fig.3a. As the maximum oscillation angle grows, the 
chaotic region in the parameter space increases on both 
sides of the second transition threshold. Then, for an¬ 
gles greater than 170°, this region approaches the first 
dynamic transition threshold where homoclinic chaos oc¬ 
curs. The absence of chaotization for maximum oscilla¬ 
tion angles less than 135° means that the system is close 
to integrable separated pendula. For larger angles the 
system is far enough from being integrable and chaos ap¬ 
pears. Such behavior is caused by interaction of dynamic 
separatrix, coinciding with LPT in the conditions of sec¬ 
ond dynamic transition, and the conventional pendulum 
separatrix. 

Summarizing, the analytical description of highly non¬ 
stationary resonance processes in a system of weakly cou¬ 
pled pendula without any restrictions on the amplitude 
of oscillations was presented for the first time. It is shown 
that such processes can be adequately described by LPTs 
corresponding to maximum possible energy exchange be¬ 
tween pendula. These regimes encircle the domains of 



FIG. 4. Poincare sections for uj — 0.65 and angles qj,max — 
37r/4 for decreasing parameter e; LPTs (red), separatrix 
(blue), a) Before first transition, e — 0.2; b) after first transi¬ 
tion, £ — 0.1; c) second transition, e m 0.0695; d) after second 
transition, £ = 0.05. 
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FIG. 5. Poincare sections for uj — 0.482 and angles qj,max — 
97t/10 for decreasing parameter e\ LPTs (red), separatrix 
(blue), a) Before first transition, e — 0.225; b) after first 
transition, £ = 0.145; c) second transition, e — 0.104; d) after 
second transition, £ = 0.03. 



(C) ‘ (d) 


FIG. 6. Time histories of the two pendula oscillations for 
<?i(0) = 37r/4 and decreasing coupling parameter e: qi(t) 
(blue), q 2 (t) (red), a) Before first transition, e = 0.2; b) 
after first transition, e = 0.1; c) second transition, e m 0.0695; 
d) after second transition, e = 0.05. 


regular motion which are determined for all initial angles 
in oscillation dynamic regime. It is also shown that man¬ 
ifestation of chaotic behavior in the considered model is 
strongly connected with a purely non-stationary dynamic 
transition. 

The obtained results can be applied in the variety of 
fields where the model of weakly coupled pendula plays a 
basic role. For example, in the application to Josephson 
junctions they correspond to a particular case in which 
both damping and external forces can be taken into ac¬ 
count in the next order approximation [3 . Therefore, the 
revealed regimes with intense inter-pendulum energy ex¬ 
change and predominant energy localization in one of the 
two pendula can be experimentally verified and exploited 
in numerous applications of Josephson junctions. 
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(C) ‘ (d) ' 

FIG. 7. Time histories of the two pendula oscillations for 
<?i(0) = 97r/10 and decreasing coupling parameter e: qi(t) 
(blue), q 2 (t) (red), a) Before first transition, e — 0.225; b) 
after first transition, e = 0.145; c) second transition, e = 
0.104; d) after second transition, e = 0.03. 
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